{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the main codes and comments on data are adapted from a notebook from Hackweek tutorial.\n",
    "The link is https://github.com/ICESAT-2HackWeek/2020_ICESat-2_Hackweek_Tutorials.git."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Libraries import\n",
    "Import all the necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import requests\n",
    "# import getpass\n",
    "# import socket\n",
    "# import json\n",
    "# import zipfile\n",
    "# import math\n",
    "# import pprint\n",
    "# import time\n",
    "\n",
    "# utility modules\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "import io\n",
    "\n",
    "# the usual suspects\n",
    "import numpy as np\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.io import shapereader\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "import matplotlib.pyplot as plt\n",
    "# import matplotlib.image as mpimg\n",
    "\n",
    "# modules you'll need if you're downloading the data\n",
    "from icepyx import icesat2data as ipd\n",
    "import shutil\n",
    "import geopandas as gpd\n",
    "\n",
    "import fiona\n",
    "import pyproj\n",
    "import h5py\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Plot ICESat-2 data\n",
    "First plot background map layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd3ed720a4dd42d59ed7c3dca652a8b6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "extent = [88.1, 89, 27.7, 28.5]\n",
    "x = [88.1443, 88.2269, 88.3643, 88.2708, 88.4037, 88.4865, 88.625, 88.5314, 88.6636, 88.7459, 88.8852, 88.7912, 88.9223, 89.0054]\n",
    "y = [27.6996, 28.5008, 27.6990, 28.5008, 27.6983, 28.5008, 27.6997, 28.5008, 27.6999, 28.5010, 27.6965, 28.5012, 27.6961, 28.5012]\n",
    "\n",
    "request = cimgt.OSM()\n",
    "\n",
    "fig = plt.figure(figsize=(9, 13))\n",
    "ax = plt.axes(projection=request.crs)\n",
    "\n",
    "ax.plot(x[0:2], y[0:2], transform=ccrs.PlateCarree(), color='red')\n",
    "ax.plot(x[2:4], y[2:4], transform=ccrs.PlateCarree(), color='red')\n",
    "ax.plot(x[4:6], y[4:6], transform=ccrs.PlateCarree(), color='red')\n",
    "ax.plot(x[6:8], y[6:8], transform=ccrs.PlateCarree(), color='red')\n",
    "ax.plot(x[8:10], y[8:10], transform=ccrs.PlateCarree(), color='red')\n",
    "ax.plot(x[10:12], y[10:12], transform=ccrs.PlateCarree(), color='red')\n",
    "ax.plot(x[12:14], y[12:14], transform=ccrs.PlateCarree(), color='red')\n",
    "\n",
    "gl = ax.gridlines(draw_labels=True, alpha=0.2)\n",
    "# gl.xlabels_top = gl.ylabels_right = False\n",
    "gl.xformatter = LONGITUDE_FORMATTER\n",
    "gl.yformatter = LATITUDE_FORMATTER\n",
    "ax.set_extent(extent)\n",
    "ax.add_image(request, 10)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c0223822fe4e46808ea1a53044ab5f40",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fiona.drvsupport.supported_drivers['kml'] = 'rw' # enable KML support which is disabled by default\n",
    "fiona.drvsupport.supported_drivers['KML'] = 'rw' # enable KML support which is disabled by default\n",
    "\n",
    "# kml data\n",
    "lake_filepath = os.getcwd()+'/sikkimlakes.kml'\n",
    "gdfl = gpd.read_file(lake_filepath) #GeoDataFrame object\n",
    "glacier_filepath = os.getcwd()+'/sikkimglaciers.kml'\n",
    "gdfg = gpd.read_file(glacier_filepath) #GeoDataFrame object\n",
    "border_filepath = os.getcwd()+'/borders.kml'\n",
    "gdfb = gpd.read_file(border_filepath)\n",
    "\n",
    "# trajectory data  x lon        y lat\n",
    "x = [88.1443, 88.2269, 88.3643, 88.2708, 88.4037, 88.4865, 88.625, 88.5314, 88.6636, 88.7459, 88.8852, 88.7912, 88.9223, 89.0054]\n",
    "y = [27.6996, 28.5008, 27.6990, 28.5008, 27.6983, 28.5008, 27.6997, 28.5008, 27.6999, 28.5010, 27.6965, 28.5012, 27.6961, 28.5012]\n",
    "\n",
    "\n",
    "# Overlay lake outline\n",
    "f, ax = plt.subplots(1, figsize=(9, 13))\n",
    "plt.title('ICESat-2 trajectory in Sikkim')\n",
    "gdfb.plot(ax=ax, facecolor='white', edgecolor='black')\n",
    "gdfl.plot(ax=ax, facecolor='lightskyblue', edgecolor='gray')\n",
    "gdfg.plot(ax=ax, facecolor='blue', edgecolor='white')\n",
    "\n",
    "plt.plot(x[0:2], y[0:2], color='red', label='477')\n",
    "plt.plot(x[2:4], y[2:4], color='olivedrab', label='241')\n",
    "plt.plot(x[4:6], y[4:6], color='gold', label='919')\n",
    "plt.plot(x[6:8], y[6:8], color='mediumslateblue', label='683')\n",
    "plt.plot(x[8:10], y[8:10], color='mediumspringgreen', label='1361')\n",
    "plt.plot(x[10:12], y[10:12], color='violet', label='1125')\n",
    "plt.plot(x[12:14], y[12:14], color='moccasin', label='416')\n",
    "plt.annotate(\"Como Chamling\", (88.19, 28.4))\n",
    "plt.legend(loc='upper right',title='track ID')\n",
    "ax.set_ylim([27.7,28.5]) #adjust to frame glacier\n",
    "ax.set_xlim([88.1,89]);\n",
    "plt.xlabel('longitude')\n",
    "plt.ylabel('latitude')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# image = mpimg.imread(\"January_original.tif\")\n",
    "# plt.imshow(image)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pre read data ATL03 ATL06 and ATL13    Reference Ground Track (RGT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc\n",
    "\n",
    "# read the IS2 data with Tyler's ATL03 reader:\n",
    "ATL03_file='ATL03_20190128163325_04770206_003_01.h5'\n",
    "IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams =read_HDF5_ATL03(ATL03_file)\n",
    "# add x_atc to the ATL03 data structure (this function adds to the LS2_ATL03_mds dictionary)\n",
    "get_ATL03_x_atc(IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7278fe5341d6443ea181f7e44242383c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "beam='gt2l'\n",
    "rgt=\"0477\"\n",
    "cycle=\"02\"\n",
    "# select the 2l beam from ATL03\n",
    "D3 = IS2_atl03_mds[beam]\n",
    "\n",
    "# create scatter plot of photon data (e.g., photon elevation vs x_atc)\n",
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(D3['heights']['x_atc'], D3['heights']['h_ph'],'k.',markersize=0.25, label='all photons')\n",
    "LMH=D3['heights']['signal_conf_ph'][:,3] >= 2\n",
    "ax.plot(D3['heights']['x_atc'][LMH], D3['heights']['h_ph'][LMH],'g.',markersize=0.5, label='flagged photons')\n",
    "h_leg=ax.legend()\n",
    "\n",
    "ax.set_xlabel('x_atc, m')\n",
    "ax.set_ylabel('h, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(type(D3['heights']))\n",
    "# print(D3['heights'].keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from readers.read_HDF5_ATL13 import read_HDF5_ATL13\n",
    "\n",
    "# # read the IS2 data with ATL13 reader:\n",
    "# ATL13_file='ATL13_20190128160153_04770201_003_01.h5'\n",
    "# IS2_atl13_mds, IS2_atl13_attrs, IS2_atl13_beams =read_HDF5_ATL13(ATL13_file)\n",
    "# D13 = IS2_atl13_mds[beam]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(D13['ht_water_surf'])\n",
    "# print(D13['ht_ortho'])\n",
    "# print(D13['ht_water_surf']-D13['ht_ortho'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ATL03 can show where the surface is, but it doesn't explicitly report a surface height. It's also a bulky product that can be slow to load. ATL06 reduces elevation data to a 20-meter posting, and gives one elevation per 20 meters (instead of hundreds in ATL03).\n",
    "To help work with ATL06 data, we'll use a simple piece of code that reads ATL06 data one beam at a time and stores it in a dictionary. The code has a default set of fields to be read from an ATL06 file, and stores the data from each field in the output dictionary. It also takes care of removing bad data, by setting values that are marked as invalid in the file to NaN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use this to read the ATL06 data that correspond to the ATL03 data in the previous cell, and will plot the ATL06 elevations on top of the ATL03 photons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the matching ATL06 file\n",
    "ATL06_file='ATL06_20190128163325_04770206_003_01.h5'\n",
    "D6=atl06_to_dict(ATL06_file,'/gt2l', index=None, epsg=3031)\n",
    "\n",
    "# plot the elevations on top of the previous axes.  You should be able to scroll up to the previous plot and see the ATL06 points.\n",
    "ax.plot(D6['x_atc'], D6['h_li'],'r.', label='ATL06')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'x_atc': 83928, 'h_li': 83928}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{field:D6[field].size for field in ['x_atc','h_li']}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ed1902aa0b0b452daa110601c7445ac7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#select the 2l beam from ATL03\n",
    "D3 = IS2_atl03_mds[beam]\n",
    "\n",
    "# find the matching ATL06 file\n",
    "D6=atl06_to_dict(ATL06_file, beam, index=None, epsg=3031)\n",
    "# create scatter plot of photon data (e.g., photon elevation vs x_atc)\n",
    "%matplotlib widget\n",
    "f1,ax = plt.subplots(num=1,figsize=(6,4))\n",
    "TEP=(D3['heights']['signal_conf_ph'][:,3] <-1)   \n",
    "ax.plot(D3['heights']['x_atc'][TEP], D3['heights']['h_ph'][TEP],'b.',markersize=0.25, label='TEP')\n",
    "BG=(D3['heights']['signal_conf_ph'][:,3] ==0)   |  (D3['heights']['signal_conf_ph'][:,3] ==1)\n",
    "ax.plot(D3['heights']['x_atc'][BG], D3['heights']['h_ph'][BG],'k.',markersize=0.25, label='Background')\n",
    "LMH=D3['heights']['signal_conf_ph'][:,3] >= 2\n",
    "ax.plot(D3['heights']['x_atc'][LMH], D3['heights']['h_ph'][LMH],'g.',markersize=0.5, label='flagged photons')\n",
    "ax.plot(D6['x_atc'], D6['h_li'],'r.', label='ATL06')\n",
    "h_leg=ax.legend()\n",
    "\n",
    "plt.title(beam+':'+os.path.basename(ATL03_file))\n",
    "\n",
    "ax.set_xlabel('x_atc, m')\n",
    "ax.set_ylabel('h, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The blue stripe is signals flagged with flag which indicates that nobody should believe that this is the surface. That is what we call the transmitter echo path arrival and these photons were fired into an optical fiber onboard the spacecraft. The photons rattled around this optical fiber for a while and came back to the detectors. \n",
    "TEP = Transmitter Echo Pulse\n",
    "Now the window of photons around the surface is much more scattered, there are some places where ATL06 is reporting heights for something that's probably not the surface, and there's a stripe of extra photons below the surface that doesn't make any obvious sense.  The new (blue) stripe of photons is the TEP or Transmitter Echo Pulse, which records a packet of photons that are looped from the transmit to the receive side of ATLAS to help monitor its impulse response.  ATL06 ignores these photons, and we should too.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gt2l\n"
     ]
    }
   ],
   "source": [
    "print(D6['beam'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from readers.read_HDF5_ATL13 import read_HDF5_ATL13\n",
    "\n",
    "# read the IS2 data with ATL13 reader:\n",
    "ATL13_file='ATL13_20190128160153_04770201_003_01.h5'\n",
    "IS2_atl13_mds, IS2_atl13_attrs, IS2_atl13_beams =read_HDF5_ATL13(ATL13_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now deal with ATL06 and ATL13 to produce glacier thicknesses and lake surface heights.\n",
    "Start of the data process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams =read_HDF5_ATL03(ATL03_file)\n",
    "D6=atl06_to_dict(ATL06_file, beam, index=None, epsg=3031)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
